# Replication package for 
# "Latent Territorial Threat and Democratic Regime Reversals"
# Johannes Karreth, Jaroslav Tir, and Douglas M. Gibler
# Forthcoming in the Journal of Peace Research, 2021
# jkarreth@ursinus.edu

# Please read "0_readme.html" before using this code.

# This file is `3_analyze_reversals.R`
# It generates the main results in the article and most of the supporting information.

rm(list = ls())

####################################
### Merge all datasets
####################################

set.seed(123)
library("tidyverse")

# Set working directory at the level of the replication folder

load(file = "Source/distpol.RData")
distpol.sm <- select(distpol.sm, -polity2)
load(file = "Source/neighbors.RData")
load(file = "Source/polity.RData")
polity.dat <- polity.dat[polity.dat$ccode != 769, ]
load(file = "Source/econ.RData")
hsigo.ipol <- read.csv("Source/HSIGOs/hsigo_cow30.csv_1945-2017.csv")
mid4_cy <- read.csv("Source/mid4_3.csv")
rivalries <- read.csv("Source/contiguity_rivalries_cy.csv")

# Fix Germany int polity.dat to not lose the pre-1990 observations
polity.dat[polity.dat$ccode == 260 & polity.dat$year >= 1955 & polity.dat$year <= 1990, ]$ccode <- 255

# Becauses the most entries are in econ.dat, start with that file

dat1 <- left_join(x = econ.dat, y = distpol.sm, by = c("ccode" = "ccode", "year" = "year"))
dat2 <- left_join(x = dat1, y = polity.dat[polity.dat$year >= 1946, ], by = c("ccode" = "ccode", "year" = "year"))
dat3 <- left_join(x = dat2, y = neighbors1946, by = c("ccode" = "ccode1", "year" = "year"))
dat4 <- left_join(x = dat3, y = hsigo.ipol, by = c("ccode" = "cown", "year" = "year"))
dat5 <- left_join(x = dat4, y = mid4_cy, by = c("ccode" = "ccode", "year" = "year"))
dat6 <- left_join(x = dat5, y = rivalries, by = c("ccode" = "ccode", "year" = "year"))

dat6$midSum <- ifelse(is.na(dat6$midSum) == TRUE, 0, dat6$midSum)
dat6$midOngoing <- ifelse(is.na(dat6$midOngoing) == TRUE, 0, dat6$midOngoing)
dat6$terrmidSum <- ifelse(is.na(dat6$terrmidSum) == TRUE, 0, dat6$terrmidSum)
dat6$terrmidOngoing <- ifelse(is.na(dat6$terrmidOngoing) == TRUE, 0, dat6$terrmidOngoing)
dat6$fatalmidSum <- ifelse(is.na(dat6$fatalmidSum) == TRUE, 0, dat6$fatalmidSum)
dat6$fatalmidOngoing <- ifelse(is.na(dat6$fatalmidOngoing) == TRUE, 0, dat6$fatalmidOngoing)
dat6$fatalTerrMidOngoing <- ifelse(dat6$fatalmidOngoing == 1 & dat6$terrmidOngoing == 1, 1, 0)

dat6$rivalry.sum <- ifelse(is.na(dat6$rivalry.sum) == TRUE, 0, dat6$rivalry.sum)
dat6$cont.rivalry.sum <- ifelse(is.na(dat6$cont.rivalry.sum) == TRUE, 0, dat6$cont.rivalry.sum)
dat6$noncont.rivalry.sum <- ifelse(is.na(dat6$noncont.rivalry.sum) == TRUE, 0, dat6$noncont.rivalry.sum)

# Territorial threat: estimated models using DR.TTpred.R, then saved as TTfm_mcmc.RData and TTpred.csv
tt <- read.csv(file  = "Work/TTpred.csv")

# Fix Germany to not lose the pre-1990 observations
tt[tt$ccode == 260 & tt$year >= 1955 & tt$year <= 1990, ]$ccode <- 255

# Create rolling mean (previous 2 & 5 years) for TT

# moving mean for that year and previous year (e.g. 5 represents the mean of that year and the four previous years)
tt2 = tt %>%
  group_by(ccode) %>%
  arrange(ccode, year) %>%
  mutate(max_terrthreatfm_mcpp_cur2 = zoo::rollmean(x = max_terrthreatfm_mcpp, 2, align = "right", fill = NA),
         max_terrthreatfm_mc1946pp_cur2 = zoo::rollmean(x = max_terrthreatfm_mc1946pp, 2, align = "right", fill = NA),
         max_terrthreatfm_mcpp_cur5 = zoo::rollmean(x = max_terrthreatfm_mcpp, 5, align = "right", fill = NA),
         max_terrthreatfm_mc1946pp_cur5 = zoo::rollmean(x = max_terrthreatfm_mc1946pp, 5, align = "right", fill = NA))
head(tt2, 75)

# moving mean for the previous years not including the current year (e.g. 5 represents the mean of the 5 previous years)
tt2 = tt2 %>%
  mutate(max_terrthreatfm_mcpp_lag1 = lag(max_terrthreatfm_mcpp, n = 1),
         max_terrthreatfm_mc1946pp_lag1 = lag(max_terrthreatfm_mc1946pp, n = 1)) %>%
  mutate(max_terrthreatfm_mcpp_pre2 = zoo::rollapply(data = max_terrthreatfm_mcpp_lag1, 
                                     width = 2, 
                                     FUN = mean, 
                                     align = "right", 
                                     fill = NA, 
                                     na.rm = TRUE),
         max_terrthreatfm_mc1946pp_pre2 = zoo::rollapply(data = max_terrthreatfm_mc1946pp_lag1, 
                                                     width = 2, 
                                                     FUN = mean, 
                                                     align = "right", 
                                                     fill = NA, 
                                                     na.rm = TRUE),
         max_terrthreatfm_mcpp_pre5 = zoo::rollapply(data = max_terrthreatfm_mcpp_lag1, 
                                                width = 5, 
                                                FUN = mean, 
                                                align = "right", 
                                                fill = NA, 
                                                na.rm = TRUE),
         max_terrthreatfm_mc1946pp_pre5 = zoo::rollapply(data = max_terrthreatfm_mc1946pp_lag1, 
                                                     width = 5, 
                                                     FUN = mean, 
                                                     align = "right", 
                                                     fill = NA, 
                                                     na.rm = TRUE))

head(tt2, 75)

# for those obs where no 2-year rolling average could be calculated, use the current value instead
tt2$max_terrthreatfm_mcpp_pre2 <- ifelse(is.na(tt2$max_terrthreatfm_mcpp_pre2) == TRUE, tt2$max_terrthreatfm_mcpp, tt2$max_terrthreatfm_mcpp_pre2)
tt2$max_terrthreatfm_mc1946pp_pre2 <- ifelse(is.na(tt2$max_terrthreatfm_mc1946pp_pre2) == TRUE, tt2$max_terrthreatfm_mc1946pp, tt2$max_terrthreatfm_mc1946pp_pre2)
tt2$max_terrthreatfm_mcpp_pre5 <- ifelse(is.na(tt2$max_terrthreatfm_mcpp_pre5) == TRUE, tt2$max_terrthreatfm_mcpp, tt2$max_terrthreatfm_mcpp_pre5)
tt2$max_terrthreatfm_mc1946pp_pre5 <- ifelse(is.na(tt2$max_terrthreatfm_mc1946pp_pre5) == TRUE, tt2$max_terrthreatfm_mc1946pp, tt2$max_terrthreatfm_mc1946pp_pre5)

dat7 <- left_join(x = dat6, y = tt2[tt2$year >= 1946, c("ccode", "year","max_terrthreatfm_mcpp", "max_terrthreatfm_mcppSD", "max_terrthreatfm_mcpp10", "max_terrthreatfm_mcpp90", "max_terrthreatfm_mcpp_cur2", "max_terrthreatfm_mcpp_cur5", "max_terrthreatfm_mcpp_lag1", "max_terrthreatfm_mcpp_pre2", "max_terrthreatfm_mcpp_pre5", "max_terrthreatfm_mc1946pp", "max_terrthreatfm_mc1946ppSD", "max_terrthreatfm_mc1946pp10", "max_terrthreatfm_mc1946pp90", "max_terrthreatfm_mc1946pp_cur2", "max_terrthreatfm_mc1946pp_cur5", "max_terrthreatfm_mc1946pp_lag1", "max_terrthreatfm_mc1946pp_pre2", "max_terrthreatfm_mc1946pp_pre5"), ], by = c("ccode" = "ccode", "year" = "year"), all.x = TRUE, all.y = FALSE)

# Remove the double Pakistan
# dat7 <- dat7[dat7$ccode != 769, ]

# A few more controls from QoG
# qog_full <- read_csv(file = "Source/QoG/qog_std_ts_jan19.csv")
qog <- read.csv(file = "Source/QoG/qog_std_ts_jan19_reduced.csv")
# qog <- qog[, c("ccodecow", "year", "ht_colonial", "ucdp_type3", "ucdp_type4")]
qog$formercolony <- ifelse(qog$ht_colonial > 0, 1, 0)
qog$formercolony <- ifelse(is.na(qog$formercolony) == TRUE, 0, qog$formercolony)
qog$frenchcolony <- ifelse(qog$formercolony == 1 & qog$ht_colonial == 6, 1, 0)
qog$britcolony <- ifelse(qog$formercolony == 1 & qog$ht_colonial == 5, 1, 0)

dat <- left_join(x = dat7, y = qog, by = c("ccode" = "ccodecow", "year" = "year"))
qog <- NULL

####################################
### 10. Lags and transformations
####################################

dat <- arrange(dat, ccode, year)
dat <- group_by(dat, ccode)
dat <- mutate(dat, aut6tran.reg.lag = lag(aut6tran.reg, n = 1), aut7tran.reg.lag = lag(aut7tran.reg, n = 1), polity2.lag = lag(polity2, n = 1), max_terrthreatfm_mcpp_ch1 = max_terrthreatfm_mcpp - lag(max_terrthreatfm_mcpp, n = 1), max_terrthreatfm_mcpp_ch2 = max_terrthreatfm_mcpp - lag(max_terrthreatfm_mcpp, n = 2), max_terrthreatfm_mcpp.lag = lag(max_terrthreatfm_mcpp, n = 1), max_terrthreatfm_mc1946pp_ch1 = max_terrthreatfm_mc1946pp - lag(max_terrthreatfm_mc1946pp, n = 1), max_terrthreatfm_mc1946pp_ch2 = max_terrthreatfm_mc1946pp - lag(max_terrthreatfm_mc1946pp, n = 2), max_terrthreatfm_mc1946pp.lag = lag(max_terrthreatfm_mc1946pp, n = 1), midPrevYr = lag(midOngoing, n = 1), fatalmidPrevYr = lag(fatalmidOngoing, n = 1), terrmidPrevYr = lag(terrmidOngoing, n = 1), fatalTerrMidPrevYr = lag(fatalTerrMidOngoing, n = 1))

dat$coldwar <- ifelse(dat$year > 1990, 1, 0)

dat$civilwar <- ifelse(dat$ucdp_type3 >= 1 | dat$ucdp_type4 >= 1, 1, 0)
dat$civilwar <- ifelse(is.na(dat$civilwar == TRUE), 0, dat$civilwar)

# Transformations of TT estimate (based on 1816-2016 sample)

# Deciles - from the whole population (not just the reversal sample)
dat$max_terrthreatfm_mcpp_q <- cut(dat$max_terrthreatfm_mcpp, breaks = quantile(dat$max_terrthreatfm_mcpp, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q <- as.numeric(dat$max_terrthreatfm_mcpp_q)
dat$max_terrthreatfm_mcpp_q_pre2 <- cut(dat$max_terrthreatfm_mcpp_pre2, breaks = quantile(dat$max_terrthreatfm_mcpp_pre2, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q_pre2 <- as.numeric(dat$max_terrthreatfm_mcpp_q_pre2)
dat$max_terrthreatfm_mcpp_q_pre5 <- cut(dat$max_terrthreatfm_mcpp_pre5, breaks = quantile(dat$max_terrthreatfm_mcpp_pre5, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q_pre5 <- as.numeric(dat$max_terrthreatfm_mcpp_q_pre5)
dat$max_terrthreatfm_mcpp.lag_q <- cut(dat$max_terrthreatfm_mcpp.lag, breaks = quantile(dat$max_terrthreatfm_mcpp.lag, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp.lag_q <- as.numeric(dat$max_terrthreatfm_mcpp.lag_q)
dat$max_terrthreatfm_mcpp_q5 <- cut(dat$max_terrthreatfm_mcpp, breaks = quantile(dat$max_terrthreatfm_mcpp, probs = seq(0, 1, by = 0.05), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q5 <- as.numeric(dat$max_terrthreatfm_mcpp_q5)

# Logarithms

dat$max_terrthreatfm_mcpp_ln <- log(dat$max_terrthreatfm_mcpp + 0.00025)
dat$max_terrthreatfm_mcppSD_ln <- log(dat$max_terrthreatfm_mcppSD + 0.00025)
dat$max_terrthreatfm_mcpp10_ln <- log(dat$max_terrthreatfm_mcpp10 + 0.00025)
dat$max_terrthreatfm_mcpp90_ln <- log(dat$max_terrthreatfm_mcpp90 + 0.00025)
dat$max_terrthreatfm_mcpp_pre2_ln <- log(dat$max_terrthreatfm_mcpp_pre2 + 0.00025)
dat$max_terrthreatfm_mcpp_pre5_ln <- log(dat$max_terrthreatfm_mcpp_pre5 + 0.00025)

# based on 1946-2016 sample

# Deciles - from the whole population (not just the reversal sample)
dat$max_terrthreatfm_mc1946pp_q <- cut(dat$max_terrthreatfm_mc1946pp, breaks = quantile(dat$max_terrthreatfm_mc1946pp, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mc1946pp_q <- as.numeric(dat$max_terrthreatfm_mc1946pp_q)
dat$max_terrthreatfm_mc1946pp_q_pre2 <- cut(dat$max_terrthreatfm_mc1946pp_pre2, breaks = quantile(dat$max_terrthreatfm_mc1946pp_pre2, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mc1946pp_q_pre2 <- as.numeric(dat$max_terrthreatfm_mc1946pp_q_pre2)
dat$max_terrthreatfm_mc1946pp_q_pre5 <- cut(dat$max_terrthreatfm_mc1946pp_pre5, breaks = quantile(dat$max_terrthreatfm_mc1946pp_pre5, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mc1946pp_q_pre5 <- as.numeric(dat$max_terrthreatfm_mc1946pp_q_pre5)
dat$max_terrthreatfm_mc1946pp.lag_q <- cut(dat$max_terrthreatfm_mc1946pp.lag, breaks = quantile(dat$max_terrthreatfm_mc1946pp.lag, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mc1946pp.lag_q <- as.numeric(dat$max_terrthreatfm_mc1946pp.lag_q)
dat$max_terrthreatfm_mc1946pp_q5 <- cut(dat$max_terrthreatfm_mc1946pp, breaks = quantile(dat$max_terrthreatfm_mc1946pp, probs = seq(0, 1, by = 0.05), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mc1946pp_q5 <- as.numeric(dat$max_terrthreatfm_mc1946pp_q5)

# Logarithms

dat$max_terrthreatfm_mc1946pp_ln <- log(dat$max_terrthreatfm_mc1946pp + 0.00025)
dat$max_terrthreatfm_mc1946ppSD_ln <- log(dat$max_terrthreatfm_mc1946ppSD + 0.00025)
dat$max_terrthreatfm_mc1946pp10_ln <- log(dat$max_terrthreatfm_mc1946pp10 + 0.00025)
dat$max_terrthreatfm_mc1946pp90_ln <- log(dat$max_terrthreatfm_mc1946pp90 + 0.00025)
dat$max_terrthreatfm_mc1946pp_pre2_ln <- log(dat$max_terrthreatfm_mc1946pp_pre2 + 0.00025)
dat$max_terrthreatfm_mc1946pp_pre5_ln <- log(dat$max_terrthreatfm_mc1946pp_pre5 + 0.00025)

# Create a standard deviation around the LOG of territorial threat
# first, create "z-scores" based on the original measure (0-100)
dat$max_terrthreatfm_mc1946pp_Z <- ifelse(dat$max_terrthreatfm_mc1946pp / dat$max_terrthreatfm_mc1946ppSD < 10, dat$max_terrthreatfm_mc1946pp / dat$max_terrthreatfm_mc1946ppSD, 10)

# then, use this z-score to calculate a SD around the logged TT score:
# log_sd = log/sd
dat$max_terrthreatfm_mc1946pp_pre2_ln_SD <- ifelse(dat$max_terrthreatfm_mc1946pp_pre2_ln / dat$max_terrthreatfm_mc1946pp_Z * (-1) > 17, 17, dat$max_terrthreatfm_mc1946pp_pre2_ln / dat$max_terrthreatfm_mc1946pp_Z * (-1))

# Other transformations

dat$demnb500_p <- dat$demnb500 * 100
dat$dem6prop.global_p <- dat$dem6prop.global * 100
dat$gdp.lag_ln <- log(dat$gdp.lag)
dat$yasdem6_ln <- log(dat$yasdem6)
dat$ysincedem6tran_ln <- log(dat$ysincedem6tran)
dat$yasdem6_sq <- dat$yasdem6^2
dat$ysincedem6tran_sq <- dat$ysincedem6tran^2
dat$yasdem6_cu <- dat$yasdem6^3
dat$ysincedem6tran_cu <- dat$ysincedem6tran^3

dat <- as.data.frame(dat)

# Restrict data to the period for which we have TT estimates (1946-2016)

# dat <- filter(dat, year <= 2016)

####################################
### 11. List of reversals
####################################

table(dat[is.na(dat$dem6 == TRUE & dat$year > 1945), ]$ccode)
table(dat[is.na(dat$dem7 == TRUE & dat$year > 1945), ]$ccode)

dat$country <- as.character(dat$country)
dat$country <- ifelse(is.na(dat$country) == TRUE, countrycode::countrycode(sourcevar = dat$ccode, origin = "cown", destination = "country.name"), dat$country)

# Inspect instances of backsliding:
dat[dat$aut6tran == 1 & is.na(dat$aut6tran) == FALSE, c("ccode", "year", "country", "polity2", "polity2.lag", "aut6tran", "yasdem6", "ysincedem6tran", "max_terrthreatfm_mcpp", "max_terrthreatfm_mc1946pp", "reversal6sample")]

dim(dat[dat$aut6tran == 1 & is.na(dat$aut6tran) == FALSE, c("ccode", "year", "country", "polity2", "polity2.lag", "aut6tran", "yasdem6", "ysincedem6tran", "max_terrthreatfm_mcpp", "max_terrthreatfm_mc1946pp", "reversal6sample")])

# 67 cases, we have TT estimates for  55 of them. 4 are dubious cases, so we code them separately (aut6tran_inclmin) below.
# Note that 2 cases are from 2017 and 2018, for which we don't have TT measure.

# ccode year            country polity2 polity2.lag aut6tran yasdem6 ysincedem6tran max_terrthreatfm_mcpp reversal6sample
# 338      41 1991              Haiti      -7           7        1       1              1          0.0099543154               1
# 346      41 1999              Haiti       2           7        1       6              5          0.0086956367               1
# 383      42 1963 Dominican Republic       0           8        1       1              1          0.0010047030               1
# 414      42 1994 Dominican Republic       5           6        1      17             16          0.0023600146               1
# 1354     91 1985           Honduras       5           6        1       3              3          0.0810141213               1
# 1386     91 2017           Honduras       5           7        1      31             28                    NA               1
# 1813    101 2006          Venezuela       5           6        1      48             48          0.0112858969               1
# 2091    135 1992               Peru      -3           8        1      12             12          0.1127001889               1
# 2364    155 1973              Chile      -7           6        1       9              9          0.0273951758               1
# 2440    160 1976          Argentina      -9           6        1       3              3          0.0248575376               1
# 2508    165 1971            Uruguay       3           8        1      19             19          0.0504355781               1
# 2933    220 1958             France       5          10        1      79             12          0.0096848377               1
# 3798    315 1947     Czechoslovakia       2          10        1      23              2          0.2044330057               1
# 4762    350 1949             Greece       4           8        1      66              5          0.3748696042               1
# 4849    352 1963             Cyprus       0           8        1       3              4                    NA               1
# 5231    365 2007             Russia       4           6        1       7              7          0.0416709560               1
# 5509    369 1993            Ukraine       5           6        1       2              3          0.0617906828               1
# 5530    369 2014            Ukraine       4           6        1      22             20          0.0049734982               1
# 5584    370 1995            Belarus       0           7        1       4              5          0.0337732504               1
# 5657    371 1995            Armenia       3           7        1       4              5          0.2741681202               1
# 6404    404 2012      Guinea-Bissau       1           6        1       7              7          0.0004088222               1
# 6532    420 1994             Gambia      -7           8        1      29             30          0.0003327091               1
# 6623    432 2012               Mali       0           7        1      20             20          0.0024676716               1
# 6899    436 1996              Niger      -6           8        1       4              4          0.0575781274               1
# 6912    436 2009              Niger      -3           6        1       9              5          0.0027097801               1
# 6919    436 2016              Niger       5           6        1      14              5          0.0058905863               1
# 7235    451 1967       Sierra Leone      -7           6        1       6              7          0.0040634476               1
# 7322    452 1981              Ghana      -7           6        1       2              2          0.0002949961               1
# 7526    475 1966            Nigeria      -7           7        1       6              7          0.0054650573               1
# 7544    475 1984            Nigeria      -7           7        1      11              5          0.1653875300               1
# 7964    500 1966             Uganda       0           7        1       4              5          0.0715231677               1
# 8232    516 2015            Burundi      -1           6        1      10             10          0.0032084465               1
# 8332    520 1969            Somalia      -7           7        1       9             10          0.1306024256               1
# 8798    551 1996             Zambia       1           6        1       5              5          0.0601165967               1
# 8949    553 2001             Malawi       4           6        1       7              7          0.0025481446               1
# 9137    570 1970            Lesotho      -9           9        1       4              5          0.0198685152               1
# 9165    570 1998            Lesotho       0           8        1       9              5          0.0290514448               1
# 9395    580 2009         Madagascar       0           7        1      17             17                    NA               1
# 9477    581 2018            Comoros      -3           9        1      14             14                    NA               1
# 9928    625 1958              Sudan      -7           8        1       2              3          0.0463876522               1
# 9939    625 1969              Sudan       2           7        1       6              4          0.1418059036               1
# 9959    625 1989              Sudan      -7           7        1       9              3          0.1131683973               1
# 10160   640 1971             Turkey      -2           8        1      19             10          0.1672700946               1
# 10169   640 1980             Turkey      -5           9        1      26              7          0.2802273844               1
# 10203   640 2014             Turkey       3           9        1      57             31          0.0106533419               1
# 10366   652 1958              Syria      NA           7        1       4              4          0.2891678401               1
# 11902   732 1961        Korea South      -7           8        1       1              1          0.4748025095               1
# 12191   770 1958           Pakistan      -7           8        1       2              1          0.2620582642               1
# 12210   770 1977           Pakistan      -7           8        1       6              4          0.1056723346               1
# 12232   770 1999           Pakistan      -6           7        1      17             11          0.3546329982               1
# 12280   771 1974         Bangladesh      -2           8        1       2              3          0.0304048451               1
# 12313   771 2007         Bangladesh      -6           6        1      18             16          0.0325967099               1
# 12341   775 1962    Myanmar (Burma)      -6           8        1      14             15          0.0746048400               1
# 12434   780 1982          Sri Lanka       5           6        1      34             35                    NA               1
# 12455   780 2003          Sri Lanka       5           6        1      36              2                    NA               1
# 12461   780 2009          Sri Lanka       5           6        1      39              3                    NA               1
# 12600   790 2002              Nepal      -6           6        1       3              3          0.0019172298               1
# 12677   800 2006           Thailand      -5           9        1      14             14          0.1124401925               1
# 12685   800 2014           Thailand      -3           7        1      17              3          0.0107667036               1
# 12777   812 1960               Laos      -1           8        1       3              2          0.1120637382               1
# 12954   820 1969           Malaysia       1          10        1      12             13          0.0291951710               1
# 12999   820 2014           Malaysia       5           6        1      18              6          0.0063237661               1
# 13021   830 1963          Singapore      NA           7        1       4              5                    NA               1
# 13715   940 2000    Solomon Islands       0           8        1      22             23                    NA               1
# 13921   950 1987               Fiji      -3           9        1      17             18                    NA               1
# 13934   950 2000               Fiji       5           6        1      18              1                    NA               1
# 13940   950 2006               Fiji      -3           6        1      20              2                    NA               1

# note 67 cases (more data than in earlier version, because now we go up to 2016)
# next: take a closer look at list before proceeding, and cross-check with previous version
# Explained in footnote 9: 
# Four additional cases were not true regime reversals but instead were artifacts of the Polity IV coding schema; we exclude them from our analysis and this discussion. The cases are: Brazil in 1947 which moved from 7 to 5 due to changes in parliamentary elections; Turkey in 1954 when electoral changes moved the country from 7 to 4; Brazil in 1961 after a 1-point change in the score due to a change in executive constraints; and Ecuador that lowered 2 points in 2007 when the president and popular vote bypassed an intransigent congress. None of these cases are labeled reversals in datasets that focus on regime changes.

####################################
### 11. Data for estimation
####################################

aut6tran.dat <- select(filter(dat, !is.na(reversal6sample) & reversal6sample == 1), ccode, year, aut6tran, max_terrthreatfm_mcpp_ch1, max_terrthreatfm_mcpp_ch2, max_terrthreatfm_mcpp, max_terrthreatfm_mcppSD, max_terrthreatfm_mcpp10, max_terrthreatfm_mcpp90, max_terrthreatfm_mcpp.lag_q, max_terrthreatfm_mcpp.lag, max_terrthreatfm_mcpp_cur2, max_terrthreatfm_mcpp_cur5, max_terrthreatfm_mcpp_pre2, max_terrthreatfm_mcpp_pre5, max_terrthreatfm_mcpp_q_pre2, max_terrthreatfm_mcpp_q_pre5, max_terrthreatfm_mcpp_pre2_ln, max_terrthreatfm_mcpp_pre5_ln, max_terrthreatfm_mcpp_ln, max_terrthreatfm_mcpp_q, max_terrthreatfm_mcppSD_ln, max_terrthreatfm_mcpp10_ln, max_terrthreatfm_mcpp90_ln, max_terrthreatfm_mc1946pp_ch1, max_terrthreatfm_mc1946pp_ch2, max_terrthreatfm_mc1946pp, max_terrthreatfm_mc1946ppSD, max_terrthreatfm_mc1946pp10, max_terrthreatfm_mc1946pp90, max_terrthreatfm_mc1946pp.lag_q, max_terrthreatfm_mc1946pp.lag, max_terrthreatfm_mc1946pp_cur2, max_terrthreatfm_mc1946pp_cur5, max_terrthreatfm_mc1946pp_pre2, max_terrthreatfm_mc1946pp_pre5, max_terrthreatfm_mc1946pp_q_pre2, max_terrthreatfm_mc1946pp_q_pre5, max_terrthreatfm_mc1946pp_pre2_ln, max_terrthreatfm_mc1946pp_pre2_ln_SD, max_terrthreatfm_mc1946pp_pre5_ln, max_terrthreatfm_mc1946pp_ln, max_terrthreatfm_mc1946pp_q, max_terrthreatfm_mc1946ppSD_ln, max_terrthreatfm_mc1946pp10_ln, max_terrthreatfm_mc1946pp90_ln, demnb500, dem6.allnb, borders, aut6tran.reg.lag, dem6prop.global, coldwar, gdppc.lag, gdppc.lag_ln, polity2.lag, aut6tran.sum, yasdem6, rcode, al_ethnic, h_j, ysincedem6tran, demnb500_p, dem6prop.global_p, gdp.lag_ln, yasdem6_ln, ysincedem6tran_ln, yasdem6_sq, ysincedem6tran_sq, yasdem6_cu, ysincedem6tran_cu, hsigo, hsigo_lag, msigo, msigo_lag, lsigo, lsigo_lag, aut6tran_inclmin, midOngoing, midSum, midPrevYr, fatalmidOngoing, fatalmidPrevYr, terrmidOngoing, terrmidPrevYr, fatalTerrMidOngoing, fatalTerrMidPrevYr, ucdp_type3, ucdp_type4, civilwar, formercolony, frenchcolony, britcolony, rivalry.sum, cont.rivalry.sum, noncont.rivalry.sum)

### Deciles - from the reversal sample only

aut6tran.dat$max_terrthreatfm_mcpp_qrs <- cut(aut6tran.dat$max_terrthreatfm_mcpp, breaks = quantile(aut6tran.dat$max_terrthreatfm_mcpp, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
aut6tran.dat$max_terrthreatfm_mcpp_qrs <- as.numeric(aut6tran.dat$max_terrthreatfm_mcpp_q)
aut6tran.dat$max_terrthreatfm_mcpp_qrs5 <- cut(aut6tran.dat$max_terrthreatfm_mcpp, breaks = quantile(aut6tran.dat$max_terrthreatfm_mcpp, probs = seq(0, 1, by = 0.05), na.rm = TRUE), include.lowest = TRUE)
aut6tran.dat$max_terrthreatfm_mcpp_qrs5 <- as.numeric(aut6tran.dat$max_terrthreatfm_mcpp_qrs5)

aut6tran.dat$max_terrthreatfm_mc1946pp_qrs <- cut(aut6tran.dat$max_terrthreatfm_mc1946pp, breaks = quantile(aut6tran.dat$max_terrthreatfm_mc1946pp, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
aut6tran.dat$max_terrthreatfm_mc1946pp_qrs <- as.numeric(aut6tran.dat$max_terrthreatfm_mc1946pp_q)
aut6tran.dat$max_terrthreatfm_mc1946pp_qrs5 <- cut(aut6tran.dat$max_terrthreatfm_mc1946pp, breaks = quantile(aut6tran.dat$max_terrthreatfm_mc1946pp, probs = seq(0, 1, by = 0.05), na.rm = TRUE), include.lowest = TRUE)
aut6tran.dat$max_terrthreatfm_mc1946pp_qrs5 <- as.numeric(aut6tran.dat$max_terrthreatfm_mc1946pp_qrs5)

# Add GDP for Guinea-Bissau in 2012 if needed
# Use value from 2011

# aut6tran.dat[aut6tran.dat$ccode == 404 & aut6tran.dat$year == 2012, ]$gdp.lag <- aut6tran.dat[aut6tran.dat$ccode == 404 & aut6tran.dat$year == 2011, ]$gdp.lag

# aut6tran.dat[aut6tran.dat$ccode == 404 & aut6tran.dat$year == 2012, ]$gdp.lag_ln <- aut6tran.dat[aut6tran.dat$ccode == 404 & aut6tran.dat$year == 2011, ]$gdp.lag_ln

# To check: Table of reversals and all predictors of the main models

select(filter(aut6tran.dat, aut6tran == 1 & !is.na(aut6tran) & !is.na(max_terrthreatfm_mcpp)), ccode, year, aut6tran, aut6tran_inclmin, max_terrthreatfm_mcpp, demnb500_p, dem6.allnb, borders, aut6tran.reg.lag, dem6prop.global_p, coldwar, civilwar, gdp.lag_ln, polity2.lag, aut6tran.sum, yasdem6_ln, al_ethnic, h_j, hsigo, max_terrthreatfm_mcpp_pre2)

# Note missing values for 
# al_ethnic (220 in 1958 [NA in QoG], 315 in 1947 [NA in QoG], 770 in 1958 [NA in QoG])

# Fill in al_ethnic from Fearon & Laitin
aut6tran.dat[aut6tran.dat$ccode == 220 & aut6tran.dat$year == 1958, ]$al_ethnic <- .272074
aut6tran.dat[aut6tran.dat$ccode == 315 & aut6tran.dat$year == 1947, ]$al_ethnic <- .504643
aut6tran.dat[aut6tran.dat$ccode == 770 & aut6tran.dat$year == 1958, ]$al_ethnic <- .5321
# aut6tran.dat[aut6tran.dat$ccode == 260, ]$al_ethnic <- .0259799


# h_j (220 in 1958 [NA in QoG], 315 in 1947 [NA in QoG], 350 in 1949 [NA in QoG], 625 in 1958 [NA in QoG], 640 in 1954 [NA in QoG], 652 in 1958 [NA in QoG], 770 in 1958 [NA in QoG])
# Can't fill in since no data pre-1960

# Repeat reversals need to be corrected for Pakistan b/c of country code issues
# aut6tran.dat[aut6tran.dat$ccode == 770 & aut6tran.dat$year == 1977, ]$aut6tran.sum <- 1
# aut6tran.dat[aut6tran.dat$ccode == 770 & aut6tran.dat$year == 1999, ]$aut6tran.sum <- 2
# aut6tran.dat[aut6tran.dat$ccode == 770 & aut6tran.dat$year == 1977, ]$aut6tran.sum_inclmin <- 1
# aut6tran.dat[aut6tran.dat$ccode == 770 & aut6tran.dat$year == 1999, ]$aut6tran.sum_inclmin <- 2


# Repeat table to check

select(filter(aut6tran.dat, aut6tran == 1 & !is.na(aut6tran) & !is.na(max_terrthreatfm_mcpp)), ccode, year, aut6tran, aut6tran_inclmin, max_terrthreatfm_mcpp, demnb500_p, dem6.allnb, borders, aut6tran.reg.lag, dem6prop.global_p, coldwar, civilwar, gdp.lag_ln, polity2.lag, aut6tran.sum, yasdem6_ln, al_ethnic, h_j, hsigo, max_terrthreatfm_mcpp_pre2)

# Problems: none

####################################
### Describe TT
####################################

TTfm.dat <- aut6tran.dat[, c("ccode", "year", "max_terrthreatfm_mcpp", "max_terrthreatfm_mcpp10", "max_terrthreatfm_mcpp90")]
# TTfmo.dat <- aut6tran.dat[, c("ccode", "year", "max_terrthreatfmo_mcpp", "max_terrthreatfmo_mcpp10", "max_terrthreatfmo_mcpp90")]
TTfm <- TTfm.dat[, c(3:5)]
# TTfmo <- TTfmo.dat[, c(3:5)]
library(reshape)
library(ggplot2)
TTfm <- melt(TTfm)
# TTfmo <- melt(TTfmo)
TTfm.p <- ggplot(data = TTfm, aes(x = value * 100, fill = variable)) + geom_histogram(alpha = 0.5, position = "identity") + theme_bw() + theme(legend.position=c(1, 0.1), legend.justification=c(1, 0.1)) + labs(fill = "Predicted territorial threat") + scale_fill_discrete(labels=c("Mean", "10th percentile", "90th percentile")) + xlab("Predicted territorial threat") + ylab("Country-years")

# ggsave(TTfm.p, file = "/Users/johanneskarreth/Documents/Dropbox/Uni/1 - Papers/36 - Democratic reversals/Data/Graphs/TTfm_hist.pdf")

# TT over time

tt_ts <- ggplot(data = dat, aes(x = year, y = max_terrthreatfm_mcpp, group = ccode, alpha = factor(reversal6sample))) + geom_line() + scale_alpha_manual(values = c(0.25, 1)) + guides(alpha = FALSE) + xlab("Year") + ylab("Territorial threat (in %)") + theme_minimal()
tt_ts <- ggplot(data = aut6tran.dat, aes(x = year, y = max_terrthreatfm_mcpp * 100, group = ccode)) + 
  geom_line(alpha = 0.25) + 
  guides(color = FALSE) + xlab("Year") + ylab("Latent territorial threat (in %)") + theme_bw() + 
  scale_x_continuous(breaks = c(1946, 1960, 1970, 1980, 1990, 2000, 2010, 2016)) + 
  geom_point(data = aut6tran.dat[aut6tran.dat$terrmidOngoing == 1, ], aes(x = year, y = max_terrthreatfm_mcpp * 100), color = "red", size = 1) + 
  annotate("text", x = 1960, y = -2, label = "Ongoing militarized interstate dispute", hjust = 0) + 
  annotate("point", x = 1959, y = -2, color = "red") + 
  annotate("text", x = 1990, y = -2, label = "Latent territorial threat estimate", hjust = 0) + 
  annotate("segment", x = 1987, xend = 1989, y = -2, yend = -2, color = "black", alpha = 0.25)

ggsave(tt_ts, file = "Output/Figures/tt_MIDs_ts.pdf", width = 10, height = 4)

# Correlation

cor(aut6tran.dat$max_terrthreatfm_mcpp, aut6tran.dat$terrmidOngoing, use = "complete.obs")
psych::biserial(x = aut6tran.dat$max_terrthreatfm_mcpp, y = aut6tran.dat$terrmidOngoing)
t.test(aut6tran.dat$max_terrthreatfm_mcpp ~ aut6tran.dat$midOngoing)

####################################
### List of reversals
####################################

# Note that gdp.lag_ln reduces our sample size and removes 6 (!) reversals
# went back to econ.dat and filled in, now improved

reversals <- model.frame(aut6tran ~ # max_terrthreatfm_mcpp + max_terrthreatfm_mcpp10 + max_terrthreatfm_mcpp90 + 
                           max_terrthreatfm_mc1946pp + max_terrthreatfm_mc1946pp10 + max_terrthreatfm_mc1946pp90 + 
                           hsigo + demnb500 + aut6tran.reg.lag + dem6prop.global + coldwar + 
                           gdppc.lag_ln + polity2.lag + aut6tran.sum + yasdem6_ln + al_ethnic + noncont.rivalry.sum + yasdem6 +
                           ccode + year, data = filter(aut6tran.dat, year <= 2016))
table(reversals$aut6tran)

# reversals <- model.frame(aut6tran ~ max_terrthreatfm_mcpp + max_terrthreatfm_mcpp10 + max_terrthreatfm_mcpp90 + demnb500_p + dem6.allnb + aut6tran.reg.lag + dem6prop.global_p + coldwar + gdp.lag_ln + polity2.lag + aut6tran.sum + yasdem6_ln + yasdem6 + ccode + year, data = aut6tran.dat)
reversals <- as.data.frame(reversals)

# Multiply TT score by 100 for percentage scores

# reversals$max_terrthreatfm_mcpp <- reversals$max_terrthreatfm_mcpp * 100
# reversals$max_terrthreatfm_mcpp10 <- reversals$max_terrthreatfm_mcpp10 * 100
# reversals$max_terrthreatfm_mcpp90 <- reversals$max_terrthreatfm_mcpp90 * 100

reversals$max_terrthreatfm_mc1946pp <- reversals$max_terrthreatfm_mc1946pp * 100
reversals$max_terrthreatfm_mc1946pp10 <- reversals$max_terrthreatfm_mc1946pp10 * 100
reversals$max_terrthreatfm_mc1946pp90 <- reversals$max_terrthreatfm_mc1946pp90 * 100

library("countrycode")
reversals$cname <- countrycode(sourcevar = reversals$ccode, origin = "cown", destination = "country.name")
reversals$cname2 <- stringi::stri_trans_totitle(reversals$cname)
reversals <- arrange(reversals, year)

library("xtable")
reversals[reversals$cname2 == "Venezuela, Bolivarian Republic Of", ]$cname2 <- "Venezuela"
reversals[reversals$cname2 == "Lao People's Democratic Republic", ]$cname2 <- "Laos"
reversals[reversals$cname2 == "Russian Federation", ]$cname2 <- "Russia"
reversals[reversals$cname2 == "Korea, Republic Of", ]$cname2 <- "South Korea"
reversals[reversals$cname2 == "Guinea-bissau", ]$cname2 <- "Guinea-Bissau"
reversals[reversals$cname2 == "Myanmar (burma)", ]$cname2 <- "Myanmar (Burma)"
rlist.tab <- reversals[reversals$aut6tran == 1 & is.na(reversals$aut6tran) == FALSE & is.na(reversals$cname2) == FALSE, c("cname2", "year", "max_terrthreatfm_mc1946pp", "polity2.lag")]
rlist.tab$year <- as.character(rlist.tab$year)
rlist.tab$polity2.lag <- as.character(rlist.tab$polity2.lag)
rlist.tab$max_terrthreatfm_mc1946pp <- round(rlist.tab$max_terrthreatfm_mc1946pp, digits = 1)
rlist <- xtable(rlist.tab, caption = "Democratic reversals, 1946-2016.", label = "tab:reversal_list", align = "llccc", digits = 1)
print(rlist, file = "Output/Tables/reversals_list.tex", caption.placement = "top", include.rownames = FALSE, booktabs = TRUE)

# Also make this into a dotplot

reversals$cyear <- paste(reversals$cname2, reversals$year)
revp.dat <- reversals[reversals$aut6tran == 1, ]

rp <- ggplot(dat = revp.dat, aes(y = max_terrthreatfm_mc1946pp, x = reorder(cyear, max_terrthreatfm_mc1946pp))) + 
  geom_hline(yintercept = quantile(aut6tran.dat$max_terrthreatfm_mc1946pp, probs = c(0.5), na.rm = TRUE) * 100, color = "darkgray") + 
  geom_pointrange(aes(y = max_terrthreatfm_mc1946pp, x = reorder(cyear, max_terrthreatfm_mc1946pp), ymin = max_terrthreatfm_mc1946pp10, ymax = max_terrthreatfm_mc1946pp90), size = 0.25)  + 
  coord_flip() + ylab("Latent territorial threat (in %)") + xlab("") + theme_bw() + 
  annotate("text", x = 2, y = 2, label = "Median territorial threat \nacross all democracy-years, 1946-2016", size = 4, colour = "darkgray", hjust = 0)

pdf("Output/Figures/tt_reversals_dotplot.pdf", width = 10, height = 10)
rp
dev.off()

rp_log <- ggplot(dat = revp.dat, aes(y = max_terrthreatfm_mc1946pp, x = reorder(cyear, max_terrthreatfm_mc1946pp))) + 
  geom_hline(yintercept = quantile(aut6tran.dat$max_terrthreatfm_mc1946pp, probs = c(0.5), na.rm = TRUE) * 100, color = "darkgray") + 
  # geom_pointrange(aes(y = max_terrthreatfm_mc1946pp, x = reorder(cyear, max_terrthreatfm_mc1946pp), ymin = max_terrthreatfm_mc1946pp10, ymax = max_terrthreatfm_mc1946pp90), size = 0.25)  + 
  geom_segment(aes(x = reorder(cyear, max_terrthreatfm_mc1946pp), xend = reorder(cyear, max_terrthreatfm_mc1946pp), y = max_terrthreatfm_mc1946pp10, yend = max_terrthreatfm_mc1946pp90), size = 1)  + 
  geom_point(aes(x = reorder(cyear, max_terrthreatfm_mc1946pp), y = max_terrthreatfm_mc1946pp), size = 2, color = "black", fill = "white", shape = 21) + 
  geom_text(dat = filter(revp.dat, max_terrthreatfm_mc1946pp < quantile(aut6tran.dat$max_terrthreatfm_mc1946pp, probs = c(0.5), na.rm = TRUE) * 100), aes(x = reorder(cyear, max_terrthreatfm_mc1946pp), y = max_terrthreatfm_mc1946pp90, label = cyear), hjust = -0.1, size = 3) + 
  geom_text(dat = filter(revp.dat, max_terrthreatfm_mc1946pp >= quantile(aut6tran.dat$max_terrthreatfm_mc1946pp, probs = c(0.5), na.rm = TRUE) * 100), aes(x = reorder(cyear, max_terrthreatfm_mc1946pp), y = max_terrthreatfm_mc1946pp10, label = cyear), hjust = 1.1, size = 3) + 
  coord_flip() + ylab("Latent territorial threat (in %, axis in logged scale)") + xlab("") + 
  theme_minimal() + 
  annotate("text", x = 2, y = 0.6, label = "Median territorial threat \nacross all democracy-years, 1946-2016", size = 4, colour = "darkgray", hjust = 0) + scale_y_log10() + scale_x_discrete(breaks = FALSE)

pdf("Output/Figures/figure1.pdf", width = 10, height = 10)
rp_log
dev.off()

ggsave(rp_log, file = "Output/Figures/figure1.eps", width = 10, height = 10)
ggsave(rp_log, file = "Output/Figures/figure1.tif", device = "tiff", dpi = 600, width = 10, height = 10, compression = "lzw", type = "cairo")
ggsave(rp_log, file = "Output/Figures/figure1.jpg", device = "jpeg", dpi = 350, width = 10, height = 10)


####################################
### Regression results in the article
####################################

##################
### standardize predictors
##################

aut6tran.dat$hsigo_std <- as.numeric(scale(aut6tran.dat$hsigo, center = TRUE, scale = TRUE))
aut6tran.dat$demnb500_std <- as.numeric(scale(aut6tran.dat$demnb500, center = TRUE, scale = TRUE))
aut6tran.dat$dem6.allnb_std <- as.numeric(scale(aut6tran.dat$dem6.allnb, center = TRUE, scale = TRUE))
aut6tran.dat$aut6tran.reg.lag_std <- as.numeric(scale(aut6tran.dat$aut6tran.reg.lag, center = TRUE, scale = TRUE))
aut6tran.dat$dem6prop.global_std <- as.numeric(scale(aut6tran.dat$dem6prop.global, center = TRUE, scale = TRUE))
aut6tran.dat$gdppc.lag_ln_std <- as.numeric(scale(aut6tran.dat$gdppc.lag_ln, center = TRUE, scale = TRUE))
aut6tran.dat$polity2.lag_std <- as.numeric(scale(aut6tran.dat$polity2.lag, center = TRUE, scale = TRUE))
aut6tran.dat$aut6tran.sum_std <- as.numeric(scale(aut6tran.dat$aut6tran.sum, center = TRUE, scale = TRUE))
aut6tran.dat$yasdem6_ln_std <- as.numeric(scale(aut6tran.dat$yasdem6_ln, center = TRUE, scale = TRUE))
aut6tran.dat$al_ethnic_std <- as.numeric(scale(aut6tran.dat$al_ethnic, center = TRUE, scale = TRUE))
aut6tran.dat$noncont.rivalry.sum_std <- as.numeric(scale(aut6tran.dat$noncont.rivalry.sum, center = TRUE, scale = TRUE))

##################
### default specification
##################

# [territorial threat measure] + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std

##################
### first, predict reversals with MID only
##################

library("rstanarm")

# fatalTerrMidOngoing, fatalTerrMidPrevYr

main1MID.eq <- aut6tran ~ fatalTerrMidOngoing + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std
main1MID.mm <- model.matrix(main1MID.eq, data = aut6tran.dat)

# Calls to stan_glm are commented out because each model takes 
# quite some time to run.
# Instead, the code loads the saved model object.
# Un-comment the call to stan_glm to fit the model on your own.

# main1MIDmcmc <- stan_glm(main1MID.eq,
#                          family = binomial(link = "logit"),
#                          data = filter(aut6tran.dat, year <= 2016), 
#                          prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
#                          prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                          QR = TRUE,
#                          cores = parallel::detectCores(),
#                          iter = 5000,
#                          control = list(adapt_delta = 0.99))
# 
# saveRDS(main1MIDmcmc, file = "Work/main1MIDmcmc.RDS")
main1MIDmcmc <- readRDS(file = "Work/main1MIDmcmc.RDS")

main2MID.eq <- aut6tran ~ fatalTerrMidPrevYr + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std
main2MID.mm <- model.matrix(main2MID.eq, data = aut6tran.dat)

# Calls to stan_glm are commented out because each model takes 
# quite some time to run.
# Instead, the code loads the saved model object.
# Un-comment the call to stan_glm to fit the model on your own.

# main2MIDmcmc <- stan_glm(main2MID.eq,
#                          family = binomial(link = "logit"),
#                          data = filter(aut6tran.dat, year <= 2016), 
#                          prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
#                          prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                          QR = TRUE,
#                          cores = parallel::detectCores(),
#                          iter = 5000,
#                          control = list(adapt_delta = 0.99))
# 
# saveRDS(main2MIDmcmc, file = "Work/main2MIDmcmc.RDS")
main2MIDmcmc <- readRDS(file = "Work/main2MIDmcmc.RDS")

# Tables

library("BayesPostEst")

main1MID_tab <- mcmcTab(main1MIDmcmc, Pr = TRUE)
main2MID_tab <- mcmcTab(main2MIDmcmc, Pr = TRUE)

mcmcReg(list(main1MIDmcmc, main2MIDmcmc),
        custom.model.names = c("Fatal territorial MID in current year", "Fatal territorial MID in previous year"),
        coefnames = list(c("Intercept", "Territorial MID", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries"),
                         c("Intercept", "Territorial MID", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries")),
        ci = 0.90,
        hpdi = TRUE,
        caption = "Posterior estimates: Territorial MIDs and democratic reversals, 1946-2016.",
        caption.above = TRUE,
        format = "tex",
        file = "Output/Tables/tab_mainMID")

rownames(main1MID_tab) <- rownames(main2MID_tab)

main_MIDft <- cbind(main1MID_tab, main2MID_tab)

main_MIDft <- data.frame(main_MIDft)
main_MIDft <- add_row(.data = main_MIDft, Median = nrow(main1MID.mm), Median.1 = nrow(main2MID.mm))
main_MIDft <- add_row(.data = main_MIDft, `Median` = 1, `Median.1` = 2, .before = 1)

main_MIDft$Variable <- c("", "Intercept", "Territorial MID", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries", "N (country-years)")

main_MIDft <- select(main_MIDft, -Lower, -Upper, -Lower.1, -Upper.1)
main_MIDft$Variable.1 <- NA

names(main_MIDft) <- c("", "Median", "Std. dev.", "Pr(Estimate)", "", "Median", "Std. dev.", "Pr(Estimate)")

library("xtable")

main_MIDxt <- xtable(main_MIDft, caption = "Posterior estimates: Territorial MIDs and democratic reversals, 1946-2016.",
                     label = "tab_mainMID", digits = 2,
                     align = "lllcccccc")

# Multicolumn add from <http://stackoverflow.com/a/28420723>

# addtorow <- list()
# addtorow$pos <- list(-1)
# addtorow$command <- paste0(paste0(' & \\multicolumn{3}{c}{', Model[!is.na(Model)], '}', collapse=''), '\\\\')

print(main_MIDxt, file = "Output/Tables/tab_mainMID.tex",
      include.colnames = TRUE, include.rownames = FALSE,
      booktabs = TRUE, caption.placement = "top", size = "scriptsize", digits = 2)

# Coefficient density plots

x1 <- as.matrix(main1MIDmcmc)[, "fatalTerrMidOngoing"]
x1.ci <- quantile(x1, probs = c(0.05, 0.95))
x2 <- as.matrix(main2MIDmcmc)[, "fatalTerrMidPrevYr"]
x2.ci <- quantile(x2, probs = c(0.05, 0.95))
dens.dat1 <- data.frame(x = density(x1)$x, y = density(x1)$y, Model = "Territorial MID in current year", ci = 0, mean = mean(x1))
dens.dat2 <- data.frame(x = density(x2)$x, y = density(x2)$y, Model = "Territorial MID in previous year", ci = 0, mean = mean(x2))
dens.dat1$ci <- with(dens.dat1, ifelse(x > x1.ci[1] & x < x1.ci[2], 1, 0))
dens.dat2$ci <- with(dens.dat2, ifelse(x > x2.ci[1] & x < x2.ci[2], 1, 0))
dens.dat <- rbind(dens.dat1, dens.dat2)

dp2 <- ggplot(data = dens.dat, aes(x = x, y = y)) + 
  geom_ribbon(data = dens.dat[dens.dat$x > 0, ], aes(x = x, ymax = y), ymin = 0, fill = "gray", alpha = 0.5) + 
  theme_bw() + xlab("Coefficient estimates") + 
  ylab("") + scale_y_continuous(breaks = NULL) + 
  geom_line(size = 0.5) + facet_wrap( ~ Model, scales = "free") 

ann_text <- data.frame(x = c(0.75, 0.7), y = c(0.1, 0.08), Model = c("Territorial MID in current year", "Territorial MID in previous year"), lab = c("70.7% > 0", "53.6% > 0"))
ann_point <- data.frame(x = c(mean(x1), mean(x2)), y = c(0, 0), Model = c("Territorial MID in current year", "Territorial MID in previous year"))
ann_segment <- data.frame(x = c(x1.ci[1], x2.ci[1]), xend = c(x1.ci[2], x2.ci[2]), y = c(0, 0), yend = c(0, 0), Model = c("Territorial MID in current year", "Territorial MID in previous year"))
dp2 <- dp2 + geom_text(data = ann_text, aes(label = lab), size = 3) + geom_point(data = ann_point, aes(x = x, y = y), size = 4) + geom_segment(data = ann_segment, aes(x = x, xend = xend, y = y, yend = yend), size = 1) + geom_vline(xintercept = 0, linetype = "longdash", size = 0.2)

ggsave(dp2, file = "Output/Figures/main12MID_density2.pdf", width = 6, height = 3)


##################
### Main models (working toward Table 1)
##################

# Note that in the article, "Model 1" refers to m3 here, the model incorporating uncertainty
# "Model 2" in the article here is m2 here, i.e. the same

main1.eq1946 <- aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std

# Calls to stan_glm are commented out because each model takes 
# quite some time to run.
# Instead, the code loads the saved model object.
# Un-comment the call to stan_glm to fit the model on your own.

# main1stan1946 <- stan_glm(main1.eq1946, 
#                       family = binomial(link = "logit"),
#                       data = filter(aut6tran.dat, year <= 2016), 
#                       prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
#                       prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                       QR = TRUE,
#                       cores = parallel::detectCores(),
#                       iter = 5000,
#                       control = list(adapt_delta = 0.99))
# 
# saveRDS(main1stan1946, file = "Work/main1stan1946.rds")
main1stan1946 <- readRDS(file = "Work/main1stan1946.rds")

# main1stan1946_infoprior <- stan_glm(main1.eq1946, 
#                           family = binomial(link = "logit"),
#                           data = filter(aut6tran.dat, year <= 2016), 
#                           prior = cauchy(location = c(0, 0, -5, 5, -5, 0, -5, 0, 0, 5, 0, 0), scale = c(rep(2.5, 12)), autoscale = TRUE),
#                           prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                           QR = TRUE,
#                           cores = parallel::detectCores(),
#                           iter = 5000,
#                           control = list(adapt_delta = 0.99))

mcmcTab(main1stan1946, Pr = TRUE)
mcmcTab(main1stan1946_infoprior, Pr = TRUE)

# JAGS

# Calls to JAGS are commented out because each model takes 
# quite some time to run.
# Instead, the code loads the saved model object.
# Un-comment the call to JAGS to fit the model on your own.

main1.1946.list <- as.list(model.frame(main1.eq1946, data = aut6tran.dat))
main1.1946.list$N <- length(main1.1946.list$aut6tran)

# main1.1946.jagsmod <- function(){
#   for(i in 1:N){
#     aut6tran[i] ~ dbern(p[i])
#     logit(p[i]) <- b[1] + 
#       b[2] * max_terrthreatfm_mc1946pp_pre2_ln[i] +
#       b[3] * hsigo_std[i] + 
#       b[4] * demnb500_std[i] + 
#       b[5] * aut6tran.reg.lag_std[i] +
#       b[6] * dem6prop.global_std[i] +
#       b[7] * coldwar[i] +
#       b[8] * gdppc.lag_ln_std[i] +
#       b[9] * polity2.lag_std[i] +
#       b[10] * aut6tran.sum_std[i] +
#       b[11] * yasdem6_ln_std[i] +
#       b[12] * al_ethnic_std[i] +
#       b[13] * noncont.rivalry.sum_std[i]
#   }
#   # MVN prior
#   # b[1:13] ~ dmnorm(b0[1:13], B0[1:13, 1:13])
#   b[1] ~ dt(0, pow(10,-2), 1)
#   for(k in 2:13){
#     b[k] ~ dt(0, pow(2.5,-2), 1)
#   }
# }
# 
# main1.1946.list$b0 <- c(rep(0, 13))
# main1.1946.list$B0 <- diag(rep(1), 13)

saveRDS(main1.1946.list, file = "Work/m1dl_forJAGS.rds")

# run m1_tojags.R instead of the code below

# main1.1946.jags.params <- c("b")
# 
# main1.1946.jags.inits <- list(list(b = rep(0, 13)),
#                               list(b = rep(0, 13)),
#                               list(b = rep(0, 13)),
#                               list(b = rep(0, 13)))
# 
# library("dclone")
# # main1.1946.fitjags <- jags.fit(data = main1.1946.list,
# #                                  params = main1.1946.jags.params,
# #                                  model = main1.1946.jagsmod,
# #                                  inits = main1.1946.jags.inits,
# #                                  n.chains = 4,
# #                                  n.iter = 500, n.burnin = 100, n.thin = 10)
# 
# cl <- makePSOCKcluster(4)
# main1.1946.fitjags.par <- jags.parfit(cl = cl,
#                                   data = main1.1946.list,
#                                   params = main1.1946.jags.params,
#                                   model = main1.1946.jagsmod,
#                                   inits = main1.1946.jags.inits,
#                                   n.chains = 4,
#                                   n.iter = 10000, n.burnin = 2000, n.thin = 10)

main1.1946.fitjags.par <- readRDS("Work/main1jags1946.rds")

mcmcTab(main1.1946.fitjags.par, Pr = TRUE)

main2.eq1946 <- aut6tran ~ max_terrthreatfm_mc1946pp_q_pre2 + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std

# Calls to stan_glm are commented out because each model takes 
# quite some time to run.
# Instead, the code loads the saved model object.
# Un-comment the call to stan_glm to fit the model on your own.

# main2stan1946 <- stan_glm(main2.eq1946, 
#                      family = binomial(link = "logit"),
#                      data = filter(aut6tran.dat, year <= 2016), 
#                      prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
#                      prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                      QR = TRUE,
#                      cores = parallel::detectCores(),
#                      iter = 5000,
#                      control = list(adapt_delta = 0.99))

# main2stan1946_infoprior <- stan_glm(main2.eq1946, 
#                           family = binomial(link = "logit"),
#                           data = filter(aut6tran.dat, year <= 2016), 
#                           prior = cauchy(location = c(0, 0, -5, 5, -5, 0, -5, 0, 0, 5, 0, 0), scale = c(rep(2.5, 12)), autoscale = TRUE),
#                           prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                           QR = TRUE,
#                           cores = parallel::detectCores(),
#                           iter = 5000,
#                           control = list(adapt_delta = 0.99))


# saveRDS(main2stan1946, file = "Work/main2stan1946.rds")
main2stan1946 <- readRDS(file = "Work/main2stan1946.rds")

mcmcTab(main2stan1946, Pr = TRUE)
mcmcTab(main2stan1946_infoprior, Pr = TRUE)

# JAGS

main2.1946.list <- as.list(model.frame(main2.eq1946, data = filter(aut6tran.dat, year <= 2016)))
main2.1946.list$N <- length(main2.1946.list$aut6tran)

saveRDS(main2.1946.list, file = "Work/m2dl_forJAGS.rds")

# run m2_tojags.R instead of the code below

main2.1946.fitjags.par <- readRDS("Work/main2jags1946.rds")

mcmcTab(main2.1946.fitjags.par, Pr = TRUE)

# xtable version

main1info_tab <- mcmcTab(main1stan1946_infoprior, Pr = TRUE)
main2info_tab <- mcmcTab(main2stan1946_infoprior, Pr = TRUE)

rownames(main1info_tab) <- rownames(main2info_tab)

main_infoft <- cbind(main1info_tab, main2info_tab)

main_infoft <- data.frame(main_infoft)
main_infoft <- add_row(.data = main_infoft, Median = as.character(nrow(model.matrix(main1stan1946_infoprior))), Median.1 = as.character(nrow(model.matrix(main2stan1946_infoprior))))
main_infoft <- add_row(.data = main_infoft, Median = "TT logged", Median.1 = "TT deciles", .before = 1)

main_infoft$Variable <- c("", "Intercept", "Territorial threat", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries", "N (country-years)")

main_infoft$Prior <- c("", "", "Neutral", "Neutral", "Negative", "Positive", "Negative", "Neutral", "Negative", "Neutral", "Neutral", "Neutral", "Positive", "Neutral", "")

main_infoft <- select(main_infoft, -Lower, -Upper, -Lower.1, -Upper.1)
main_infoft$Variable.1 <- NA
main_infoft <- main_infoft[ , c(1, 9, 2:8)]

names(main_infoft) <- c("", "Prior", "Median", "Std. dev.", "Pr(Estimate)", "", "Median", "Std. dev.", "Pr(Estimate)")

library("xtable")

main_infoxt <- xtable(main_infoft, caption = "Posterior estimates, based on informative prior distributions: Territorial threat and democratic reversals, 1946-2016. Territorial threat estimates based on the 1946-2016 period.",
                     label = "tab_main1946_infoprior", digits = 2,
                     align = "llllcccccc")

# Multicolumn add from <http://stackoverflow.com/a/28420723>

addtorow <- list()
addtorow$pos <- list(-1)
addtorow$command <- paste0(paste0(' & \\multicolumn{3}{c}{', Model[!is.na(Model)], '}', collapse=''), '\\\\')

print(main_infoxt, file = "Output/Tables/tab_main1946_infoprior.tex",
      add.to.row = addtorow, include.colnames = TRUE, include.rownames = FALSE,
      booktabs = TRUE, caption.placement = "top", size = "scriptsize", digits = 2)

##################
### Model 3 -  TT with uncertainty
##################

m3mf <- model.frame(aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln + max_terrthreatfm_mc1946pp_pre2_ln_SD + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std + ccode + year, data = aut6tran.dat)

m3dl <- as.list(m3mf)
m3dl$N <- nrow(m3mf)
m3dl$N.ccode <- length(unique(m3dl$ccode))
m3dl$N.yr <- length(unique(m3dl$year))
m3dl$yearID <- as.numeric(as.ordered(m3dl$year))
m3dl$ccodeID <- as.numeric(as.ordered(m3dl$ccode))
# m3dl$max_terrthreatfm_mc1946ppSD <- ifelse(m3dl$max_terrthreatfm_mc1946ppSD == 0, min( m3dl$max_terrthreatfm_mc1946ppSD[m3dl$max_terrthreatfm_mc1946ppSD != min(m3dl$max_terrthreatfm_mc1946ppSD)] ), m3dl$max_terrthreatfm_mc1946ppSD)

saveRDS(m3dl, file = "Work/m3dl_forJAGS.rds")

# Using JAGS

# run m3_tojags.R instead of the code below

main3.1946.fitjags.par <- readRDS(file = "Work/main3jags1946.rds")

# Table

m3_tab <- mcmcTab(main3.1946.fitjags.par, ci = 0.9, Pr = TRUE)[, c(1, 2, 3, 6)]
m1_tab <- mcmcTab(main1.1946.fitjags.par, ci = 0.9, Pr = TRUE)[, c(1, 2, 3, 6)]
m2_tab <- mcmcTab(main2.1946.fitjags.par, ci = 0.9, Pr = TRUE)[, c(1, 2, 3, 6)]

m3_ft <- cbind(m3_tab, m1_tab, m2_tab)
Model <- c("TT logged (with uncertainty)", rep(NA, 2), "TT logged", rep(NA, 2), "TT deciles", rep(NA, 2))
N <- c(as.character(length(m3dl$aut6tran)), rep(NA, 2), as.character(length(m3dl$aut6tran)), rep(NA, 2), as.character(length(m3dl$aut6tran)), rep(NA, 2))
rownames(m3_ft) <- m3_tab[, 1]
m3_ft <- data.frame(m3_ft)
m3_ft <- rbind(m3_ft, N)
m3_ft$Variable <- c("Intercept", "Territorial threat", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries", "Country-years")

m3_ft <- m3_ft[c(2:13, 1, 14), c(1, 2, 3, 4, 6, 7, 8, 10, 11, 12)]
rownames(m3_ft) <- NULL
names(m3_ft) <- c("", "Mean", "Std. dev.", "p-direction", "Mean", "Std. dev.", "p-direction", "Mean", "Std. dev.", "p-direction")

library("xtable")

m3_xt <- xtable(m3_ft, caption = "Posterior estimates: Territorial threat (including measurement uncertainty) and democratic reversals, 1946-2016.",
                label = "tab_main1946_uncertainty", digits = 2,
                align = "llccccccccc")

# Multicolumn add from <http://stackoverflow.com/a/28420723>

addtorow <- list()
addtorow$pos <- list(-1)
addtorow$command <- paste0(paste0(' & \\multicolumn{3}{c}{', Model[!is.na(Model)], '}', collapse=''), '\\\\')

print(m3_xt, file = "Output/Tables/tab_main1946_uncertainty_312.tex", 
      add.to.row = addtorow, include.colnames = TRUE, include.rownames = FALSE,
      booktabs = TRUE, caption.placement = "top", size = "footnotesize", digits = 2)

### Plot posterior density

m3mcmc <- as.matrix(main3.1946.fitjags.par)
m1mcmc <- as.matrix(main1.1946.fitjags.par)
m2mcmc <- as.matrix(main2.1946.fitjags.par)

x1 <- m3mcmc[, 2]
x1.ci <- quantile(x1, probs = c(0.05, 0.95))
x2 <- m1mcmc[, 2]
x2.ci <- quantile(x2, probs = c(0.05, 0.95))
x3 <- m2mcmc[, 2]
x3.ci <- quantile(x3, probs = c(0.05, 0.95))
dens.dat1 <- data.frame(x = density(x1)$x, y = density(x1)$y, Model = "(1) TT logged (with uncertainty)", ci = 0, mean = mean(x1))
dens.dat2 <- data.frame(x = density(x2)$x, y = density(x2)$y, Model = "(3) TT logged", ci = 0, mean = mean(x2))
dens.dat3 <- data.frame(x = density(x3)$x, y = density(x3)$y, Model = "(2) TT deciles", ci = 0, mean = mean(x3))
dens.dat1$ci <- with(dens.dat1, ifelse(x > x1.ci[1] & x < x1.ci[2], 1, 0))
dens.dat2$ci <- with(dens.dat2, ifelse(x > x2.ci[1] & x < x2.ci[2], 1, 0))
dens.dat3$ci <- with(dens.dat3, ifelse(x > x3.ci[1] & x < x3.ci[2], 1, 0))
# dens.dat <- rbind(dens.dat1, dens.dat2, dens.dat3)
dens.dat <- rbind(dens.dat1, dens.dat3)

# dens.dat$Model_f <- factor(dens.dat$Model, levels = c("TT logged (with uncertainty)", "TT deciles"))

dp2 <- ggplot(data = dens.dat, aes(x = x, y = y)) + 
  geom_ribbon(data = dens.dat[dens.dat$x > 0, ], aes(x = x, ymax = y), ymin = 0, fill = "lightgray") + 
  theme_bw() + xlab("Estimated logit coefficient for territorial threat variable") + ylab("") + 
  scale_y_continuous(breaks = NULL) + geom_line(size = 0.5) + 
  facet_wrap( ~ (Model), scales = "free") 

ann_text <- data.frame(x = c(0.19,0.13), y = c(0.5, 0.5), Model = c("(1) TT logged (with uncertainty)", "(2) TT deciles"), lab = c("95.2% > 0", "95.5% > 0"))
ann_point <- data.frame(x = c(mean(x1), mean(x3)), y = c(0, 0), Model = c("(1) TT logged (with uncertainty)", "(2) TT deciles"))
ann_segment <- data.frame(x = c(x1.ci[1], x3.ci[1]), xend = c(x1.ci[2], x3.ci[2]), y = c(0, 0), yend = c(0, 0), Model = c("(1) TT logged (with uncertainty)", "(2) TT deciles"))
dp2 <- dp2 + geom_text(data = ann_text, aes(label = lab), size = 4) + geom_point(data = ann_point, aes(x = x, y = y), size = 4) + geom_segment(data = ann_segment, aes(x = x, xend = xend, y = y, yend = yend), size = 1) + geom_vline(xintercept = 0, linetype = "longdash", size = 0.2)

dp2 <- dp2 + annotate(geom = "text", x = 0.2, y = 4, label = "Pr(Estimate)", hjust = 0) + 
  annotate(
    geom = "segment", x = 0.25, y = 3.75, xend = 0.15, yend = 0.8, 
    # curvature = .3, 
    arrow = arrow(length = unit(2, "mm"))
  )

ggsave(dp2, file = "Output/Figures/figure2.pdf", width = 6, height = 3)
ggsave(dp2, file = "Output/Figures/figure2.eps", width = 6, height = 3)
ggsave(dp2, file = "Output/Figures/figure2.tif", device = "tiff", dpi = 600, , width = 6, height = 3, compression = "lzw", type = "cairo")
ggsave(dp2, file = "Output/Figures/figure2.jpg", device = "jpeg", dpi = 350, , width = 6, height = 3)


##################
### Predicted probabilities
##################

# Model 1 (TT logged)

# threat_range <- seq(from = min(model.matrix(main1stan1946)[, 2]), to = max(model.matrix(main1stan1946)[, 2]), length.out = 20)
# 
# main1stan1946_ppobs <- mcmcObsProb(modelmatrix = model.matrix(main1stan1946), 
#                                 mcmcout = as.matrix(main1stan1946), 
#                                 xcol = 2, 
#                                 xrange = threat_range, 
#                                 link = "logit", ci = c(0.05, 0.95), fullsims = FALSE)

main1_1946_mm <- model.matrix(main1.eq1946, data = aut6tran.dat)

threat_range1 <- seq(from = min(main1_1946_mm[, 2]), to = max(main1_1946_mm[, 2]), length.out = 20)

main1jags1946_ppobs <- mcmcObsProb(modelmatrix = main1_1946_mm, 
                                   mcmcout = as.matrix(main1.1946.fitjags.par)[sample(1:nrow(as.matrix(main1.1946.fitjags.par)), size = 1000), ], 
                                   xcol = 2, 
                                   xrange = threat_range1, 
                                   link = "logit", ci = c(0.05, 0.95), fullsims = FALSE)

pp1.plot <- ggplot(data = main1jags1946_ppobs, aes(x = x, y = median_pp * 100))
pp1.plot <- pp1.plot + geom_line() + 
  geom_line(aes(x = x, y = lower_pp * 100), linetype = "dashed") + 
  geom_line(aes(x = x, y = upper_pp * 100), linetype = "dashed")
pp1.plot <- pp1.plot + theme(axis.text.x = element_blank()) + xlab("Territorial Threat (%)") + 
  ylab("Probability of democratic reversal (in %)") + theme_bw() + 
  scale_x_continuous(limits = c(min(main1_1946_mm[, 2]), to = max(main1_1946_mm[, 2])), breaks = seq(min(main1_1946_mm[, 2]), to = max(main1_1946_mm[, 2]), length.out = 8), labels = round(c(exp(seq((min(threat_range1)), (max(threat_range1)), length.out = 8)) * 100), digits = 1))

main1_1946_tt <- data.frame(tt = (main1_1946_mm[, 2]))

pp1.hist <- ggplot(data = main1_1946_tt, aes(x = (tt))) + geom_histogram() + 
  theme_bw() + xlab("Territorial Threat (%)") + ylab("Obs.") + 
  scale_x_continuous(limits = c(min(main1_1946_mm[, 2]), to = max(main1_1946_mm[, 2])), breaks = seq(min(main1_1946_mm[, 2]), to = max(main1_1946_mm[, 2]), length.out = 8), labels = round(c(exp(seq((min(threat_range1)), (max(threat_range1)), length.out = 8)) * 100), digits = 1))

library("cowplot")

pdf("Output/Figures/main1_pp.pdf", width = 4, height = 8)
plot_grid(pp1.plot, pp1.hist, align = "v", ncol = 1, rel_heights = c(4, 1))
dev.off()

# Model 2 (TT deciles)

main2_1946_mm <- model.matrix(main2.eq1946, data = aut6tran.dat)

threat_range2 <- sort(unique(main2_1946_mm[, 2]))

main2jags1946_ppobs <- mcmcObsProb(modelmatrix = main2_1946_mm, 
                                   mcmcout = as.matrix(main2.1946.fitjags.par)[sample(1:nrow(as.matrix(main2.1946.fitjags.par)), size = 1000), ], 
                                   xcol = 2, 
                                   xrange = threat_range2, 
                                   link = "logit", ci = c(0.05, 0.95), fullsims = FALSE)

pp2.plot <- ggplot(data = main2jags1946_ppobs, aes(x = factor(x), y = median_pp * 100)) + 
  # geom_boxplot(notch = TRUE, outlier.shape = 1) + 
  geom_point() + 
  geom_segment(aes(xend = factor(x), y = lower_pp * 100, yend = upper_pp * 100)) + 
  xlab("Territorial Threat (Deciles)") + ylab("Probability of democratic reversal (in %)") + 
  # stat_summary(fun.y = "mean", geom = "point", shape = 20, size = 3, color = "red") + 
  theme_bw()

main2_1946_tt <- data.frame(tt = (main2_1946_mm[, 2]))

pp2.hist <- ggplot(data = main2_1946_tt, aes(x = (factor(tt)))) + geom_histogram(stat = "count") + 
  theme_bw() + xlab("Territorial Threat (Deciles)") + ylab("Obs.")

pdf("Output/Figures/main2_pp.pdf", width = 4, height = 8)
plot_grid(pp2.plot, pp2.hist, align = "v", ncol = 1, rel_heights = c(4, 1))
dev.off()

# all four in one

pdf("Output/Figures/main12_pp.pdf", width = 9, height = 6)
plot_grid(pp1.plot, pp2.plot, pp1.hist, pp2.hist, align = "v", ncol = 2, rel_heights = c(4, 1))
dev.off()


#########################
###  First differences
#########################

fd1 <- mcmcFD(modelmatrix = main1_1946_mm, 
       mcmcout = as.matrix(main1.1946.fitjags.par), 
       link = "logit", ci = c(0.05, 0.95), 
       percentiles = c(0.1, 0.9), 
       fullsims = TRUE)

fd2 <- mcmcFD(modelmatrix = main2_1946_mm, 
              mcmcout = as.matrix(main2.1946.fitjags.par), 
              link = "logit", ci = c(0.05, 0.95), 
              percentiles = c(0.1, 0.9), 
              fullsims = TRUE)

tt_fd1 <- fd1[, 1]
tt_fd1.ci <- quantile(tt_fd1, probs = c(0.05, 0.95))
tt_fd2 <- fd2[, 1]
tt_fd2.ci <- quantile(tt_fd2, probs = c(0.05, 0.95))

dens.dat1 <- data.frame(x = density(tt_fd1)$x, y = density(tt_fd1)$y, Model = "TT logged", ci = 0, mean =  mean(tt_fd1))
dens.dat2 <- data.frame(x = density(tt_fd2)$x, y = density(tt_fd2)$y, Model = "TT deciles", ci = 0, mean = mean(tt_fd2))

dens.dat1$ci <- with(dens.dat1, ifelse(x > tt_fd1.ci[1] & x < tt_fd1.ci[2], 1, 0))
dens.dat2$ci <- with(dens.dat2, ifelse(x > tt_fd2.ci[1] & x < tt_fd2.ci[2], 1, 0))

dens.dat <- rbind(dens.dat1, dens.dat2)

dp2 <- ggplot(data = dens.dat, aes(x = x, y = y)) + 
  geom_ribbon(data = dens.dat[dens.dat$x > 0, ], aes(x = x, ymax = y), ymin = 0, fill = "gray", alpha = 0.5) + 
  theme_bw() + xlab("Change in Pr(Democratic reversal) if TT moves from 10th to 90th percentile") + ylab("") + 
  scale_y_continuous(breaks = NULL) + geom_line(size = 0.5) + 
  facet_wrap( ~ Model, scales = "free") + xlim(-0.005, 0.02)

ann_text <- data.frame(x = c(0.01, 0.01), y = c(50, 50), Model = c("TT logged", "TT deciles"), lab = c("95.2% > 0", "95.5% > 0"))
ann_point <- data.frame(x = c(mean(tt_fd1), mean(tt_fd2)), y = c(0, 0), Model = c("TT logged", "TT deciles"))
ann_segment <- data.frame(x = c(tt_fd1.ci[1], tt_fd2.ci[1]), xend = c(tt_fd1.ci[2], tt_fd2.ci[2]), y = c(0, 0), yend = c(0, 0), Model = c("TT logged", "TT deciles"))

mult_format <- function() {
  function(x) format(100*x,digits = 2) 
}

dp2 <- dp2 + geom_text(data = ann_text, aes(label = lab), size = 4, hjust = 0) + 
  geom_point(data = ann_point, aes(x = x, y = y), size = 4) + 
  geom_segment(data = ann_segment, aes(x = x, xend = xend, y = y, yend = yend), size = 1) + 
  geom_vline(xintercept = 0, linetype = "longdash", size = 0.2) + 
  scale_x_continuous(labels = mult_format())

ggsave(dp2, file = "Output/Figures/main12_fd.pdf", width = 6, height = 3)


##################
### Model fit
##################

# Separation plot

library("separationplot")

main1.fit <- as.numeric(plogis(main1_1946_mm %*% apply(as.matrix(main1.1946.fitjags.par), MARGIN = 2, FUN = mean)))
  
main2.fit <- as.numeric(plogis(main2_1946_mm %*% apply(as.matrix(main2.1946.fitjags.par), MARGIN = 2, FUN = mean)))

separationplot(pred = main1.fit, actual = model.frame(main1.eq1946, data = aut6tran.dat)[, 1], type = "line", line = TRUE, shuffle = TRUE, show.expected = FALSE, zerosfirst = TRUE, BW = FALSE, file = "Output/Figures/main1_sep.pdf", width = 3, height = 1)
separationplot(pred = main2.fit, actual = model.frame(main2.eq1946, data = aut6tran.dat)[, 1], type = "line", line = TRUE, shuffle = TRUE, show.expected = FALSE, zerosfirst = TRUE, BW = FALSE, file = "Output/Figures/main2_sep.pdf", width = 3, height = 1)

main1_curves <- mcmcRocPrcGen(modelmatrix = main1_1946_mm, 
              mcmcout = as.matrix(main1.1946.fitjags.par)[sample(1:nrow(as.matrix(main1.1946.fitjags.par)), size = 1000), ], 
              modelframe = model.frame(main1.eq1946, data = aut6tran.dat), 
              curves = TRUE,
              link = "logit", fullsims = FALSE)

main1_roc <- ggplot(data = main1_curves$roc_dat, aes(x = x, y = y)) +
  geom_line() + 
  geom_abline(intercept = 0, slope = 1, color = "gray") + 
  labs(title = "ROC curve for Model 2 (TT logged)") + 
  xlab("1 - Specificity") + 
  ylab("Sensitivity") + theme_bw() + 
  annotate(geom = "text", x = 0.75, y = 0.25, label = "Area under curve: 0.86")

main2_curves <- mcmcRocPrcGen(modelmatrix = main2_1946_mm, 
                              mcmcout = as.matrix(main2.1946.fitjags.par)[sample(1:nrow(as.matrix(main2.1946.fitjags.par)), size = 1000), ], 
                              modelframe = model.frame(main2.eq1946, data = aut6tran.dat), 
                              curves = TRUE,
                              link = "logit", fullsims = FALSE)

main2_roc <- ggplot(data = main2_curves$roc_dat, aes(x = x, y = y)) +
  geom_line() + 
  geom_abline(intercept = 0, slope = 1, color = "gray") + 
  labs(title = "ROC curve for Model 3 (TT deciles)") + 
  xlab("1 - Specificity") + 
  ylab("Sensitivity") + theme_bw() + 
  annotate(geom = "text", x = 0.75, y = 0.25, label = "Area under curve: 0.87")

pdf("Output/Figures/main12_roc.pdf", width = 9, height = 4.5)
grid.arrange(main1_roc, main2_roc, ncol = 2, nrow = 1)
dev.off()

####################################
### Matching
####################################

library("Matching")

match.eq <- aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln + max_terrthreatfm_mc1946pp_q_pre2 + hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std

match_y <- model.frame(match.eq, data = filter(aut6tran.dat, year >= 1946))[, "aut6tran"]
match_Xmat <- model.frame(match.eq, data = aut6tran.dat)[, 4:14]

aut6tran.match <- Match(Tr = match_y, X = match_Xmat, M = 40)
aut6tran.mb <- MatchBalance(match_y ~ as.matrix(match_Xmat), match.out = aut6tran.match)

# Balance plots

library("cobalt")
aut6tran.baltab <- bal.tab(aut6tran.match, aut6tran ~ hsigo_std + demnb500_std + aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + noncont.rivalry.sum_std, data = model.frame(match.eq, data = aut6tran.dat))

v <- data.frame(old = c("hsigo_std", "demnb500_std", "aut6tran.reg.lag_std", "dem6prop.global_std", "coldwar", "gdppc.lag_ln_std", "polity2.lag_std", "aut6tran.sum_std", "yasdem6_ln_std", "al_ethnic_std", "noncont.rivalry.sum_std"),
                new = c("HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries"))

balplot <- love.plot(aut6tran.baltab,
                     # stats = c("mean.diffs", "statistics"),
                     labels = TRUE,
                     grid = TRUE,
                     colors = c("#386cb0", "#000000"),
                     shapes = c("triangle", "circle"),
                     stars = "std",
                     var.names = v) + xlab("Mean differences") + 
  labs(title = "Covariate balance before and after matching", caption = "* indicates standardized variables")

ggsave(balplot, file = "Output/Figures/balplot.pdf", width = 6, height = 3)


matched.obs <- c(aut6tran.match$index.treated, aut6tran.match$index.control)
# matched.obs <- c(unique(aut6tran.match$matches[, 1]), unique(aut6tran.match$matches[, 2]))

aut6tran.matched.dat <- model.frame(match.eq, data = aut6tran.dat)
aut6tran.matched.dat$index <- c(1:nrow(aut6tran.matched.dat))

aut6tran.matched.dat <- aut6tran.matched.dat[aut6tran.matched.dat$index %in% matched.obs, ]

matched.mod1 <- glm(aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln, data = aut6tran.matched.dat, family = "binomial", x = TRUE)
summary(matched.mod1)

matched.mod2 <- glm(aut6tran ~ max_terrthreatfm_mc1946pp_q_pre2, data = aut6tran.matched.dat, family = "binomial", x = TRUE)
summary(matched.mod2)

detach("package:Matching")
detach("package:MASS")

library("rstanarm")

# Calls to stan_glm are commented out because each model takes 
# quite some time to run.
# Instead, the code loads the saved model object.
# Un-comment the call to stan_glm to fit the model on your own.

# main1.matched.mcmc <- stan_glm(aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln, 
#                                data = aut6tran.matched.dat, 
#                                family = binomial(link = "logit"),
#                                prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
#                                prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                                cores = parallel::detectCores(),
#                                iter = 5000,
#                                control = list(adapt_delta = 0.99))
# 
# main2.matched.mcmc <- stan_glm(aut6tran ~ max_terrthreatfm_mc1946pp_q_pre2, 
#                                data = aut6tran.matched.dat, 
#                                family = binomial(link = "logit"),
#                                prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
#                                prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
#                                cores = parallel::detectCores(),
#                                iter = 5000,
#                                control = list(adapt_delta = 0.99))

library("BayesPostEst")
mcmcTab(main2.matched.mcmc, Pr = TRUE)

# mcmcReg(list(main1.matched.mcmc, main2.matched.mcmc),
#         custom.model.names = c("TT logged, matched sample", "TT deciles, matched sample"),
#         coefnames = list(c("Intercept", "Territorial threat (logged)"),
#                          c("Intercept", "Territorial threat (deciles)")),
#         ci = 0.90,
#         hpdi = FALSE,
#         caption = "Posterior estimates, based on matched samples: Territorial threat and democratic reversals, 1946-2016. Territorial threat estimates based on the 1946-2016 period.",
#         caption.above = TRUE,
#         format = "tex",
#         file = "Output/Tables/tab_main1946_matched")

# xtable version

main1match_tab <- mcmcTab(main1.matched.mcmc, Pr = TRUE)
main2match_tab <- mcmcTab(main2.matched.mcmc, Pr = TRUE)

rownames(main2match_tab) <- rownames(main1match_tab)

mainmatchft <- cbind(main1match_tab, main2match_tab)

mainmatchft <- data.frame(mainmatchft)
mainmatchft <- add_row(.data = mainmatchft, Median = as.character(nrow(model.matrix(main1.matched.mcmc))), Median.1 = as.character(nrow(model.matrix(main2.matched.mcmc))))
mainmatchft <- add_row(.data = mainmatchft, Median = "TT logged", Median.1 = "TT deciles", .before = 1)

mainmatchft$Variable <- c("", "Intercept", "Territorial threat", "N (country-years)")

mainmatchft <- dplyr::select(mainmatchft, -Lower, -Upper, -Lower.1, -Upper.1)
mainmatchft$Variable.1 <- NA

names(mainmatchft) <- c("", "Median", "Std. dev.", "Pr(Estimate)", "", "Median", "Std. dev.", "Pr(Estimate)")

library("xtable")

mainmatchxt <- xtable(mainmatchft, caption = "Posterior estimates from matched sample: Territorial threat and democratic reversals, 1946-2016. Territorial threat estimates based on the 1946-2016 period. See Figure \ref{fig:balplot} for variables on which the sample was matched (with replacement).",
                     label = "tab_main1946_matched", digits = 2,
                     align = "lllcccccc")

# Multicolumn add from <http://stackoverflow.com/a/28420723>

addtorow <- list()
addtorow$pos <- list(-1)
addtorow$command <- paste0(paste0(' & \\multicolumn{3}{c}{', Model[!is.na(Model)], '}', collapse=''), '\\\\')

print(mainmatchxt, file = "Output/Tables/tab_main1946_matched.tex",
      add.to.row = addtorow, include.colnames = TRUE, include.rownames = FALSE,
      booktabs = TRUE, caption.placement = "top", size = "scriptsize", digits = 2)

##################
### BMA
##################

aut6tran.bma.dat <- model.frame(aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln + max_terrthreatfm_mc1946pp_q_pre2 + hsigo_std + demnb500_std + 
                                  aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + 
                                  polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + 
                                  noncont.rivalry.sum_std, data = aut6tran.dat)

names(aut6tran.bma.dat) <- c("Reversal", "Territorial threat (logged)", "Territorial threat (Deciles)", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries")

library("randomForest")
aut6tran.rf1 <- randomForest(y = factor(aut6tran.bma.dat[, 1]),
                            x = aut6tran.bma.dat[, c(2, 4:14)],
                            importance = TRUE)

importance(aut6tran.rf1)

pdf("Output/Figures/main1_varImpPlot.pdf", width = 12, height = 6)
varImpPlot(aut6tran.rf1, main = "")
dev.off()

aut6tran.rf2 <- randomForest(y = factor(aut6tran.bma.dat[, 1]),
                             x = aut6tran.bma.dat[, c(3:14)],
                             importance = TRUE)

importance(aut6tran.rf2)
varImpPlot(aut6tran.rf2)

library("BMA")
aut6tran.bma1 <- bic.glm(y = factor(aut6tran.bma.dat[, 1]),
                         x = aut6tran.bma.dat[, c(2, 4:14)], 
                         glm.family = "binomial")
summary(aut6tran.bma1)

#################################################
# Years as democracy since democratization (not overall)
#################################################

main1.longeq1946 <- aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln + hsigo_std + demnb500_std + 
  aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + 
  polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + 
  noncont.rivalry.sum_std + ysincedem6tran + yasdem6 + ccode + year
main1.mf <- model.frame(main1.longeq1946, data = aut6tran.dat)
m1.dat <- main1.mf
# m1.dat <- merge(x = main1.mf, y = aut6tran.dat[, c("ccode", "year", "ysincedem6tran_sq")], by = c("ccode", "year"), all.x = TRUE, all.y = FALSE)
# m1.dat$ysincedem6tran <- sqrt(m1.dat$ysincedem6tran_sq)
# m1.dat$yasdem6 <- exp(m1.dat$yasdem6_ln)

table(m1.dat$aut6tran)
table(m1.dat[m1.dat$aut6tran == 1, ]$ysincedem6tran)
table(m1.dat[m1.dat$aut6tran == 1, ]$yasdem6)
hist(m1.dat[m1.dat$aut6tran == 1, ]$ysincedem6tran)
quantile(m1.dat[m1.dat$aut6tran == 1, ]$ysincedem6tran, probs = seq(0, 1, by = 0.1), na.rm = TRUE)
quantile(m1.dat[m1.dat$aut6tran == 1, ]$yasdem6, probs = seq(0, 1, by = 0.1), na.rm = TRUE)

ysincedem6tran.dat <- as.data.frame(table(m1.dat[m1.dat$aut6tran == 1, ]$ysincedem6tran))
yasdem6.dat <- as.data.frame(table(m1.dat[m1.dat$aut6tran == 1, ]$yasdem6))
ysincedem6tran.dat$Var1 <- as.numeric(levels(ysincedem6tran.dat$Var1))
yasdem6.dat$Var1 <- as.numeric(levels(yasdem6.dat$Var1))
years.dat <- merge(x = ysincedem6tran.dat, y = yasdem6.dat, by = "Var1", all.x = TRUE, all.y = TRUE)
names(years.dat) <- c("ReversalCases", "Years since democratization", "Years as democracy")

years.dat <- reshape2::melt(years.dat, id.vars = "ReversalCases")

years.p <- ggplot(data = years.dat, aes(x = ReversalCases, y = value, fill = variable)) + geom_bar(stat = "identity", position = "dodge", size = 2) + scale_y_continuous(name = "Reversals") + scale_fill_manual(values=c("#fdc086", "#386cb0")) + xlab("")+ theme_bw() + xlab("Years") + theme(legend.position = c(0.85, 0.8), legend.justification = c(1, 1)) + labs(fill = "")

ggsave(filename = "Output/Figures/reversals_barchart.pdf", width = 7, height = 3)

#################################################
# Summary statistics
#################################################

library("pastecs")
library("knitr")

summary.mf <- model.frame(aut6tran ~ max_terrthreatfm_mc1946pp_pre2_ln + max_terrthreatfm_mc1946pp_q_pre2 + 
                            hsigo_std + demnb500_std + 
                            aut6tran.reg.lag_std + dem6prop.global_std + coldwar + gdppc.lag_ln_std + 
                            polity2.lag_std + aut6tran.sum_std + yasdem6_ln_std + al_ethnic_std + 
                            noncont.rivalry.sum_std, data = aut6tran.dat)
summary.dat <- t(stat.desc(summary.mf))
summaryind.dat <- summary.dat[, c( "mean", "std.dev", "min", "max", "nbr.val")]
summaryind.tab <- round(summaryind.dat, digits = 2)
rownames(summaryind.tab) <- c("Democratic reversals", "Territorial threat (logged)", "Territorial threat (Deciles)", "HSIGO memberships", "Perc. democratic within 500km", "Reversals in region", "Perc. democratic (global)", "Post-Cold War", "GDPpc (t-1, logged)", "Polity (t-1)", "Previous reversals", "Years as democracy (logged)", "Ethnic fractionalization", "Non-contiguous rivalries")
colnames(summaryind.tab) <- c("Mean/Proportion", "Std. Dev.", "Min.", "Max.", "N (country-years)")
summaryind.tab[, 5] <- as.character(summaryind.tab[, 5])
kable(summaryind.tab, format = "markdown", row.names = TRUE, caption = "Summary statistics")

summaryind.xtab <- xtable(summaryind.tab, caption = "Summary statistics for models of democratic reversals.", label = "tab_sumstat", digits = 2, align = "lccccc")
print(summaryind.xtab, file = "Output/Tables/tab_sumstat.tex", include.colnames = TRUE, include.rownames = TRUE, booktabs = TRUE, caption.placement = "top", size = "footnotesize", digits = 2)

#################################################
# End of script
#################################################

# > sessionInfo()
# R version 4.0.3 (2020-10-10)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Big Sur 10.16
# 
# Matrix products: default
# LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] dclone_2.3-0       Matrix_1.2-18      coda_0.19-4        knitr_1.30         pastecs_1.3.21     cobalt_4.2.4       Matching_4.9-7    
# [8] separationplot_1.3 foreign_0.8-80     MASS_7.3-53        Hmisc_4.4-2        Formula_1.2-4      survival_3.2-7     lattice_0.20-41   
# [15] RColorBrewer_1.1-2 cowplot_1.1.0      BayesPostEst_0.2.1 rstanarm_2.21.1    Rcpp_1.0.6         xtable_1.8-4       countrycode_1.2.0 
# [22] reshape_0.8.8      forcats_0.5.0      stringr_1.4.0      dplyr_1.0.2        purrr_0.3.4        readr_1.4.0        tidyr_1.1.2       
# [29] tibble_3.0.4       ggplot2_3.3.2      tidyverse_1.3.0   
# 
# loaded via a namespace (and not attached):
#   [1] utf8_1.1.4           tidyselect_1.1.0     lme4_1.1-26          htmlwidgets_1.5.3    grid_4.0.3           miscTools_0.6-26    
# [7] jtools_2.1.1         munsell_0.5.0        codetools_0.2-16     statmod_1.4.35       DT_0.16              miniUI_0.1.1.1      
# [13] withr_2.3.0          colorspace_2.0-0     highr_0.8            rstudioapi_0.13      ROCR_1.0-11          stats4_4.0.3        
# [19] gbRd_0.4-11          bayesplot_1.7.2      Rdpack_2.1           labeling_0.4.2       rstan_2.21.2         mnormt_2.0.2        
# [25] clarkeTest_0.1.0     farver_2.0.3         vctrs_0.3.6          generics_0.1.0       xfun_0.20            R6_2.5.0            
# [31] markdown_1.1         VGAM_1.1-4           bitops_1.0-6         assertthat_0.2.1     promises_1.1.1       scales_1.1.1        
# [37] nnet_7.3-14          texreg_1.37.5        gtable_0.3.0         processx_3.4.5       sandwich_3.0-0       rlang_0.4.10        
# [43] splines_4.0.3        checkmate_2.0.0      rjags_4-10           broom_0.7.3          inline_0.3.17        yaml_2.2.1          
# [49] reshape2_1.4.4       abind_1.4-5          modelr_0.1.8         unmarked_1.0.1       threejs_0.3.3        crosstalk_1.1.0.1   
# [55] backports_1.2.1      httpuv_1.5.4         rsconnect_0.8.16     DAMisc_1.6.2         tools_4.0.3          psych_2.0.12        
# [61] runjags_2.0.4-6      optiscale_1.2        ellipsis_0.3.1       raster_3.4-5         R2WinBUGS_2.1-21     ggridges_0.5.2      
# [67] plyr_1.8.6           base64enc_0.1-3      ps_1.5.0             prettyunits_1.1.1    rpart_4.1-15         zoo_1.8-8           
# [73] cluster_2.1.0        haven_2.3.1          fs_1.5.0             survey_4.0           magrittr_2.0.1       data.table_1.13.6   
# [79] openxlsx_4.2.3       colourpicker_1.1.0   lmtest_0.9-38        reprex_0.3.0         effects_4.2-0        tmvnsim_1.0-2       
# [85] matrixStats_0.57.0   hms_0.5.3            shinyjs_2.0.0        mime_0.9             evaluate_0.14        shinystan_2.5.0     
# [91] rio_0.5.16           jpeg_0.1-8.1         readxl_1.3.1         gridExtra_2.3        rstantools_2.1.1     compiler_4.0.3      
# [97] bdsmatrix_1.3-4      V8_3.4.0             crayon_1.3.4         minqa_1.2.4          StanHeaders_2.21.0-7 htmltools_0.5.0     
# [103] later_1.1.0.1        RcppParallel_5.0.2   lubridate_1.7.9.2    DBI_1.1.0            ggmcmc_1.5.0         dbplyr_2.0.0        
# [109] boot_1.3-25          AICcmodavg_2.3-1     car_3.0-10           cli_2.2.0            mitools_2.4          rbibutils_2.0       
# [115] insight_0.11.1       igraph_1.2.6         pkgconfig_2.0.3      sp_1.4-4             R2jags_0.6-1         xml2_1.3.2          
# [121] dygraphs_1.1.1.6     plm_2.2-5            rvest_0.3.6          snakecase_0.11.0     callr_3.5.1          digest_0.6.27       
# [127] janitor_2.0.1        rmarkdown_2.6        cellranger_1.1.0     htmlTable_2.1.0      maxLik_1.4-6         curl_4.3            
# [133] shiny_1.5.0          gtools_3.8.2         nloptr_1.2.2.2       lifecycle_0.2.0      nlme_3.1-149         jsonlite_1.7.2      
# [139] carData_3.0-4        fansi_0.4.1          pillar_1.4.7         GGally_2.0.0         loo_2.4.1            fastmap_1.0.1       
# [145] httr_1.4.2           pkgbuild_1.2.0       glue_1.4.2           xts_0.12.1           zip_2.1.1            png_0.1-7           
# [151] shinythemes_1.1.2    pander_0.6.3         stringi_1.5.3        caTools_1.18.0       latticeExtra_0.6-29 
